OpenHDF5
Allie Weill
June 20, 2016
About this code
#load libraries
library(raster)## Loading required package: splibrary(rhdf5)
library(rgdal)## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3Let’s Open a File!!
f<-"../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"
# view h5 structure
h5ls(f)## group name otype dclass
## 0 / ATCOR_Input_File H5I_DATASET STRING
## 1 / ATCOR_Processing_Log H5I_DATASET STRING
## 2 / Aerosol Optical Depth H5I_DATASET INTEGER
## 3 / Aspect H5I_DATASET FLOAT
## 4 / Cast Shadow H5I_DATASET INTEGER
## 5 / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6 / Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7 / Illumination Factor H5I_DATASET INTEGER
## 8 / Path Length H5I_DATASET FLOAT
## 9 / Processing Version H5I_DATASET STRING
## 10 / Reflectance H5I_DATASET INTEGER
## 11 / Shadow_Processing_Log H5I_DATASET STRING
## 12 / Sky View Factor H5I_DATASET INTEGER
## 13 / Skyview_Processing_Log H5I_DATASET STRING
## 14 / Slope H5I_DATASET FLOAT
## 15 / Slope_Aspect_Processing_Log H5I_DATASET STRING
## 16 / Solar Azimuth Angle H5I_DATASET FLOAT
## 17 / Solar Zenith Angle H5I_DATASET FLOAT
## 18 / Surface Elevation H5I_DATASET FLOAT
## 19 / Visibility Index Map H5I_DATASET INTEGER
## 20 / Water Vapor Column H5I_DATASET INTEGER
## 21 / coordinate system string H5I_DATASET STRING
## 22 / flightAltitude H5I_DATASET FLOAT
## 23 / flightHeading H5I_DATASET FLOAT
## 24 / flightTime H5I_DATASET FLOAT
## 25 / fwhm H5I_DATASET FLOAT
## 26 / map info H5I_DATASET STRING
## 27 / to-sensor azimuth angle H5I_DATASET FLOAT
## 28 / to-sensor zenith angle H5I_DATASET FLOAT
## 29 / wavelength H5I_DATASET FLOAT
## dim
## 0 1
## 1 1
## 2 544 x 578 x 1
## 3 544 x 578 x 1
## 4 544 x 578 x 1
## 5 544 x 578 x 1
## 6 544 x 578 x 1
## 7 544 x 578 x 1
## 8 544 x 578 x 1
## 9 1
## 10 544 x 578 x 426
## 11 1
## 12 544 x 578 x 1
## 13 1
## 14 544 x 578 x 1
## 15 1
## 16 1 x 1
## 17 1 x 1
## 18 544 x 578 x 1
## 19 544 x 578 x 1
## 20 544 x 578 x 1
## 21 1
## 22 5732053
## 23 5732053
## 24 5732053
## 25 426 x 1
## 26 1
## 27 544 x 578 x 1
## 28 544 x 578 x 1
## 29 426 x 1Import Spatial Information
# import spatial info
mapInfo <- h5read(f,
"map info",
read.attributes = TRUE)
mapInfo## [1] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## attr(,"Description")
## [1] "Basic Map information for envi style programs"Grab Reflectance Materials
## read in reflectance data attributes
reflInfo <- h5readAttributes(file = f,
name = "Reflectance")
## define scale factor
scaleFactor <- reflInfo$`Scale Factor`
## define no data value
noDataValue <- as.numeric(reflInfo$`data ignore value`)
str(noDataValue)## num 15000Import Data Dims
# Open the file for viewing
fid <- H5Fopen(f)
# open the reflectance data
did <- H5Dopen(fid, "Reflectance")
# grab the dataset dimensions (column, row, band)
sid <- H5Dget_space(did)
dims <- H5Sget_simple_extent_dims(sid)$size
dims## [1] 544 578 426# close all open connections
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)Read in Reflectance Data
# extract slce of H5 file
b56 <- h5read(f,
"Reflectance",
index = list(1:dims[1],1:dims[2],56))Convert to Matrix
# convert to matrix
b56 <- b56[,,1]Plot Data
## let's plot some data
## first change from scientific notation to decimals
options("scipen"=100, "digits"=4)
## now try an image
image(b56)
## log transform it--too dark, can't see well otherwise
image(log(b56),main="log transformed data")
## histogram
hist(b56,
col="springgreen",
main="Distribution of Reflectance Values \nBand 56")
## Data Clean Up
# assign no data value to object (already created noDataValue earlier)
# set all reflectance values = 15,000 to NA
b56[b56 == noDataValue] <- NA
# apply scale factor (created earlier)
# reflectance values between 0-1 (the valid range)
b56 <- b56/scaleFactor
# view distribution of reflectance values
hist(b56,
main="Distribution with NoData Value Considered\nData Scaled")
Transposing (flipping) the data
# Because the data import as column, row but we require row, column in R, we need to transpose x and y values in order for our final image to plot properly
b56 <- t(b56)
image (log (b56),
main="Band 56\nTransposed Values")
## Define spatial extent
# We can extract the upper left-hand corner coordinates.
# the numbers as position 4 and 5 are the UPPER LEFT CORNER (x,y)
mapInfo <- strsplit(mapInfo, ",")
mapInfo <- unlist(mapInfo)
mapInfo## [1] "UTM" "1" "1"
## [4] "325963.0" "4103482.0" "1.0000000000e+000"
## [7] "1.0000000000e+000" "11" "North"
## [10] "WGS-84" "units=Meters"# grab the X,Y left corner coordinate ensure the format is numeric with as.numeric()
xMin <- as.numeric(mapInfo[4])
yMax <- as.numeric(mapInfo[5])
# we can get the x and y resolution from this string too
res <- c(mapInfo[2],mapInfo[3])
res <- as.numeric(res)
# finally calculate the xMax value and the yMin value from the dimensions
# we grabbed above. The xMax is the left corner + number of columns* resolution.
xMax <- xMin + (dims[1]*res[1])
yMin <- yMax - (dims[2]*res[2])
# also note that x and y res are the same (1 meter) Now, define the raster extent define the extent (left, right, top, bottom)
rasExt <- extent(xMin, xMax,yMin,yMax)
# now we can create a raster and assign it its spatial extent
## CRS is a coordinate reference system
b56r <- raster(b56,
crs=CRS("+init=epsg:32611"))
# assign CRS
extent(b56r) <- rasExt
# view raster object attributes
b56r## class : RasterLayer
## dimensions : 578, 544, 314432 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4102904, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0.0033, 0.5678 (min, max)# plot the new image
plot(b56r, main="Raster for Lower Teakettle \nBand 56")
Import NEON Functions
library(devtools)
#install_github("lwasser/neon-aop-package/neonAOP")
library(neonAOP)
??read_band
??open_band
b55 <- open_band(f,
band = 55,
epsg = 32611)
b55## class : RasterLayer
## dimensions : 578, 544, 314432 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4102904, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0.0029, 0.5622 (min, max)## plot data
plot(b55, main="Raster for Lower Teakettle Band 55")
# import several bands
bands <- c(58, 34, 19)
## create a raster stack
RGBstack <- create_stack(f,
bands = bands,
epsg = 32611)
plot(RGBstack)
# plot RGB stack
plotRGB(RGBstack,
stretch="lin")
#cir image
bandsCIR <- c(90, 34, 19)
## create a raster stack
CIRstack <- create_stack(f,
bands = bandsCIR,
epsg = 32611)
plot(CIRstack)
# plot RGB stack
plotRGB(CIRstack,
stretch="lin")